// copyright:no copyright,use it if you like
// license boilerplate:BSD 3-Clause
// filename:EquationSolver.hpp
// author:Tempest(https://gitee.com/Tempest56890)
// date:2022-10-3
// brief:implement the bisection method, Newton’s method, and the secant method.

#ifndef __TEMPEST_EQUATIONSOLVER__
#define __TEMPEST_EQUATIONSOLVER__

#include <iostream>
#include <cmath>
#include <limits>
#include <stdlib.h>

const double _eps=std::numeric_limits<double>::epsilon();

// implement the funciton of BaseClass.
// provides several functions for derived classes to use.
// example:
//    (Temporarily omitted.)
class EquationSolver {
public:
    // Defines a virtual function solve().
    virtual double solve()=0;
    void exchange(double *a, double *b) {
        double c=*a;
        *a=*b;
        *b=c;
    }
    int sgn(double x) {
        if(x==0) return 0;
        else if(x<0) return -1;
        else return 1;
    }
};

// implement the funciton of BaseClass.
// implements derivative functions.
// example:
//    (Temporarily omitted.)
class Function {
public:
    // Defines a virtual function operator().
    virtual double operator()(const double &x) const =0;
    // First-Order derivative function is defined.
    virtual double First_Order_diff(const double &x) const {
        return ((*this)(x+100*_eps)-(*this)(x-100*_eps))/(_eps*200);
    }
    // Second-Order derivative function is defined.
    virtual double Second_Order_diff(const double &x) const {
        double left_diff=((*this)(x)-(*this)(x-200*_eps))/(_eps*200);
        double right_diff=((*this)(x+200*_eps)-(*this)(x))/(_eps*200);
        return (right_diff-left_diff)/(_eps*200);
    }
};

// implement the funciton of derived class Bisection.
// The function is to solve the root of the equation.
// example:
//    (Temporarily omitted.)
class Bisection:public EquationSolver {
private:
    double a,b,delta;
    Function &f;
    int M;
public:
    Bisection(Function &f,double a,double b,double delta,int M)
    : f(f),a(a),b(b),delta(delta),M(M) {}
    double solve() {
        double h,u=f(a),v=f(b),c=0.0,w=0.0;
        int k=1;
        for(;k<=M;k++) {
            h=b-a;
            c=a+h/2;
            w=f(c);
            if(fabs(h)<delta||fabs(w)<_eps) {
                break;
            } else if (sgn(w)!=sgn(u)) {
                b=c;
                v=w;
            } else {
                a=c;
                u=w;
            }
        }
        std::cout<<"Bisection Method completed."<<std::endl;
        std::cout<<"The iterations is "<<k<<"."<<" The root r = "<<c<<"."<<" Now f(r) = "<<w<<"."<<std::endl;
        return c;
    }
};

// implement the funciton of derived class Newton.
// The function is to solve the root of the equation.
// example:
//    (Temporarily omitted.)
class Newton:public EquationSolver {
private:
    double x0;
    int M;
    Function &f;
public:
    Newton(Function &f,double x0,int M)
    : f(f),x0(x0),M(M) {}
    double solve() {
        double x=x0,u=0.0;
        int k=0;
        for(;k<=M;k++) {
            u=f(x);
            if(fabs(u)<_eps) {
                break;
            }
            x=x-u/f.First_Order_diff(x);
        }
        std::cout<<"Newton Method completed."<<std::endl;
        std::cout<<"The iterations is "<<k<<"."<<"The root r = "<<x<<"."<<" Now f(r) = "<<u<<"."<<std::endl;
        return x;
    }
};

// implement the funciton of derived class Secant.
// The function is to solve the root of the equation.
// example:
//    (Temporarily omitted.)
class Secant:public EquationSolver {
private:
    double x0,x1,delta;
    Function &f;
    int M;
public:
    Secant(Function &f,double x0,double x1,double delta,int M)
    : f(f),x0(x0),x1(x1),delta(delta),M(M) {}
    double solve() {
        double x_new=x1,x_old=x0,u=f(x_new),v=f(x_old),s;
        int k=2;
        for(;k<=M;k++) {
            if(fabs(u)>fabs(v)) {
                exchange(&x_old,&x_new);
                exchange(&u,&v);
            }
            s=(x_new-x_old)/(u-v);
            x_old=x_new;
            v=u;
            x_new=x_new-u*s;
            u=f(x_new);
            if(fabs(x_new-x_old)<delta||fabs(u)<_eps) {
                break;
            }
        }
        std::cout<<"Secant Method completed."<<std::endl;
        std::cout<<"The iterations is "<<(k-2)<<"."<<" The root r = "<<x_new<<"."<<" Now f(r) = "<<u<<"."<<std::endl;
        return x_new;
    }
};

#else
#endif
